Formation and switching of defect dipoles in acceptor doped lead titanate: 
A kinetic model based on first-principles calculations 
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The formation and field-induced switching of defect dipoles in acceptor doped lead titanate is described by 
a kinetic model representing an extension of the well established Arlt-Neumann model [Ferroelectrics 76, 303 
(1987)]. Energy barriers for defect association and reorientation of oxygen vacancy-dopant (Cu and Fe) com- 
plexes are obtained from first-principles calculations and serve as input data for the kinetic coefficients in the 
rate equation model. The numerical solution of the model describes the time evolution of the oxygen vacancy 
distribution at different temperatures and dopant concentrations in the presence or absence of an alternating ex- 
ternal field. We predict the characteristic time scale for the alignment of all defect dipoles with the spontanenous 
polarization of the surrounding matrix. In this state the defect dipoles act as obstacles for domain wall motion 
and contribute to the experimentally observed aging. Under cycling conditions the fully aligned configuration is 
perturbed and a dynamic equilibrium is established with defect dipoles in parallel and anti-parallel orientation 
relative to the spontaneous polarization. This process can be related to the deaging behavior of piezoelectric 



ceramics. 



I. INTRODUCTION 

Aging phenomena, namely the gradual change of phys- 
ical properties with time, are observed in almost all 
ferroelectrics^^ In some acceptor doped barium titanate 
(BaTiOa) and lead zirconate titanate (PZT) ceramics aging 
goes along with an increasing shift of the hysteresis along the 
axis of the electrical field giving rise to an internal bias field. 14 
In the past, several plausible models have been developed to 
intepret the occurence of bias fields and aging phenomena 
in ferrocelectrics in terms of domain splitting 2 , space-charge 
formation, 4 electronic charge trapping, 7,8 , ionic drift 11 and re- 
orientation of defect dipoles^^ 7 - 

In acceptor ("hard") doped ferroelectrics transition metals 
usually substitute the B-site (Ti or Zr in PZT) and tend to 
bind strongly to oxygen vacancies. These acceptor center- 
oxygen vacancy associates form electric and elastic defect 
dipoles such as charged (Fe^ Ti %")' or (Cu£. )Ti -F ") x ^ 
which contribute to the overall polarization in a ferroelec- 
tric compound- 7 *^— and can be aligned either parallel, anti- 
parallel, or perpendicular to the polarization of the surround- 
ing material as shown schematically in Fig.Q] 

In the parelectric state, defect dipoles of different orien- 
tation are energetically equivalent, whereas they have a pre- 
ferred orientation in a polar matrix. Arlt and Neumann^ 5 , 
have attributed the occurence of internal bias fields to the 
switching of defect dipoles and described the transient orien- 
tation of dipoles by a kinetic model. As quantitative data on 
the energy landscape for these defect dipoles was unavailable 
at the time they relied on a very simple electrostatic estimate 
of the energy difference. 26 The energetic asymmetry between 
the parallel and anti-parallel dipoles obtained in this fashion 
for BaTiOa was about 30 meV and thus much smaller than the 
energy differences calculated more recently by first-principles 
methods for PbTi03, I8 ~ 20 which revealed that the energetic 
asymmetry is actually as large as the barriers for oxygen mi- 
gration. Only recently, Marton and Elsasser— showed that in 



Fe-doped lead titanate the barrier for reorientation sensitively 
depends on the position of the migrating oxygen vacancy with 
respect to the iron atom and the surrounding spontaneous fer- 
roelectric polarization. During fast field cycling, the defect- 
dipoles are expected not to change orientation, because the 
characteristic rate for oxygen jumps around the acceptor cen- 
ter should be lower than the domain switching process. Ex- 
perimentally, Zhang et al. m followed the dynamics of (Mnji- 
Vq) x dipoles in barium titanate by electron paramagnetic res- 
onance studies and found support for the so-called "defect 
symmetry principle", which assumes that non-switching de- 
fect dipoles impose a restoring force for reversible domain 
switching. 15 Jakes et al. could show that in Fe 3+ doped PZT 
defect dipoles are not preferentially located at domain walls 
but within the domains. 27 Morozov et. al. studied aging- 
deaging process in hard PZT ceramics using the harmonic 
analysis of polarization response under switching conditions 
and concluded that two or more mechanisms are responsible 
for domain stabilization. 1 1 Activation energies of about 0.6 eV 
were attributed to short-range charge hopping, which could be 
due to local reorientation of microdipoles. 

Since the switching dynamics of defect dipoles depends on 
the electric and thermal energy provided to change polariza- 
tion direction, the contribution of dipole reorientation can only 
be reliably assessed if realistic numbers for the migration and 
association energies are available, which allow to quantita- 
tively model the switching dynamics of defect dipoles in a 
comprehensive way. 

The objective of the present work is to develop a kinetic 
model that captures the formation of defect dipoles as well as 
their reorientation both in the absence and presence of elec- 
tric fields^ Cu and Fe-doped lead titanate are considered as 
representative examples and the energy landscape for oxygen 
vacancy migration in these materials is obtained using first- 
principles calculations. We consider both free oxygen vacan- 
cies and oxygen vacancies associated with Fe or Cu. Starting 
from a statistical distribution the majority of oxygen vacan- 
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cies is initially unbound. Over time vacancies are captured by 
impurity atoms and subsequently converted into the lowest en- 
ergy configuration, which corresponds to a defect dipole that 
is aligned parallel to the macroscopic polarization (Mb — V& 
in Fig. [TJ. While the exact time scales for these processes are 
dependent on dopant type, concentration, and temperature our 
results demonstrate that the ground state is reached within sec- 
onds at temperature slightly above room temperature and thus 
that already the pistine material can be considered as "aged". 
In the presence of an oscillating external field our model pre- 
dicts that a gradual reorientation of defect dipoles leads to a 
dynamic equilibrium, in which the parallel and anti-parallel 
configurations occur with equal probability. This is in accord 
with the experimental observation of deaging by the applica- 
tion of AC fields.-^ 

This paper is organized as follows. First, we describe the 
kinetic model and discuss its features. This is followed in 
Sect. |nl] by a description of the first-principles calculations 
that were carried to determine the model parameters. In 
Sect. |IV] we apply the kinetic model to study vacancy redis- 
tribution as a function of temperature and impurity concentra- 
tion both in the absence and presence of an oscillating external 
electric field. The implications of the present findings for ag- 
ing and fatigue are discussed in Sect. [V] and conclusions are 
summarized in Sect. [VI] 



II. KINETIC MODEL 

In this section we formulate a kinetic model that describes 
the redistribution of oxygen vacancies between different types 
of sites as a function of time. It captures the temperature, im- 
purity concentration and frequency dependence of this pro- 
cess within a mean-field approximation. Figure [T]provides an 
overview of the different types of oxygen vacancies that are 
taken into account by this model. 

In general, the temporal variation of the concentration of 
vacancies of type i can be described by a rate equation 
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(1) 



where the first term on the right hand side accounts for the 
"loss" of vacancies of type i while the second term describes 
the "gain" due to vacancy jumps from sites of type j to sites 
of type i. At typical device operation temperatures near 300 K 
the creation or annihilation of vacancies at surfaces or inter- 
faces is negligible and the total concentration of vacancies can 
be assumed as constant, 



(2) 



The rate at which vacancies of type i jump onto sites of type 
j is given by 



Kij = u exp - 



where Vq is the attempt frequency and AG*^- 7 is the free en- 
ergy of migration encountered by a vacancy jumping from a 
site of type i to a site of type j. The attempt frequency is for 
all jumps approximated by the frequency of the lowest optical 
mode at V, u « 2 THz^ 

The probability <&jj for a vacancy to jump from a site of 
type i to a site of type j is given by the fraction of sites of type 
j in the first nearest neighbor shell of sites of type i. Using a 
simple mapping to index different defect configurations, Mg- 
V cl ^ (1), M B -V c2 ^ (2), M B -V ab ^ (3), V c ^ (4), and V ab ^ 
(5) [see Fig. [TJg) for examples] and taking into account the 
geometry of the lattice (see Fig. [TJ the following probability 
matrix is obtained 
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where a = 6/m- Here /m is the fraction of B sites which 
have been replaced by impurity atoms. The recurrence of the 
factor eight in Eq. (01 results from the number of oxygen sites 
in the second neighbor shell of any given oxygen site, while 
the factor six stems from the number of oxygen sites in the 
first neighbor shell of a B site. Introducing W,j = $yKy- 
and Vjj = 5ij J2k Wit, Eq. (Q]i can be rewritten in a conve- 



nient matrix form 
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which in this work has been numerically 32 using an adaptive 
time step algorithm for stiff differential equations. 33 

It should be noted that the model does not take into account 
the possibility of two or more oxygen vacancies associating 
with the same impurity atom. Both experiments and calcula- 
tions indicate, however, that this is unlikely to occur for the 
dopants considered in the present work. Similarly the possi- 
bility that two or more impurity atoms form an aggregate can 
be ruled out based on experimental evidenced 

In the present form the model does not include any con- 
straints to allow for the number of oxygen vacancies to be 
larger than the number of impurity atoms or vice versa. This 
situation can, however, be implemented rather easily by solv- 
ing the kinetic model in steps. For instance, consider a case 
in which the vacancy concentration is [Vo] = 0.01 and the 
dopant/impurity concentration is [M] = 0.005. The sum of 
the relative concentrations of complexed vacancies [M]/[Vb] 
can, therefore, not exceed c max = 0.5 > c\ + ci + c-j. 
Starting from some initial distribution, one solves the kinetic 
model until c\ + C2 + C3 equals c max . At this point all 
impurities are complexed with vacancies, and one can "re- 
move" the free vacancy concentrations C4 and C5 from the 
model. This is achieved by reducing the 5x5 matrices in 
Eqs. (Q]|3]l and Eq. (@) to 3 x 3 matrices, only keeping ele- 
ments G {1, 2, 3}. The opposite scenario, in which the 
number of impurity /dopant atoms exceeds the number of free 
vacancies, can be implemented in a similar fashion. For the 
sake of clarity and because the key conclusions of this work 
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FIG. 1. (a) Unit cell of the ideal tetragonal perovskite (ABO3) lattice, (b-c) Free oxygen vacancies, (d-f) Complexed oxygen vacancies 
on nearest-neighbor sites of B-site impurities. The dark blue arrow on the left indicates the direction of the lattice polarization whereas the 
smaller yellow arrows represent the orientation of the defect dipoles. (g) Schematic representation of the possible paths for oxygen vacancy 
migration. Some migration barriers, AGJ^ J , are exemplarily indicated. In (b-g) Only the B and O sites are shown. Oxygen vacancies jump 
along the vertices of the BOe octahedron. Blue and yellow circles represent native and impurity atoms on B sites, respectively. The dark and 
light red circles indicate oxygen atoms in the first nearest neighbor shell of native and impurity atoms, respectively. 



are unaffected by these conditions, we do not consider any of 
these cases in the remainder part of this paper. 

Applying the model to a specific material requires knowl- 
edge of the energy differences between various vacancy con- 
figurations as well as migration energies. To provide these pa- 
rameters, we have carried out first-principles calculations that 
are described in the following section. At this level we ne- 
glected the vibrational entropy contribution to the migration 
barriers and approximated AG m » AE m . 



III. FIRST-PRINCIPLES CALCULATIONS 
A. Computational parameters 

The barriers for oxygen vacancy migration in pure as well 
Cu and Fe-doped lead titanate were calculated within density 
functional theory (DFT) using the Vienna ab-initio simula- 
tion package. 34 The potentials due to the ions and the core 
electrons were represented by the projector-augmented wave 
method. 35 The 5d electrons of Pb, the 3s and 3p electrons of 
Ti as well as the 3p electrons of Fe and Cu were treated as 
part of the valence. The exchange-correlation potential was 
represented using the local spin density approximation,^. Su- 
percells containing 2x2x4 unit cells equivalent to 80 atoms 
were employed and the Brillouin zone was sampled using a 
2x2x2 Monkhorst-Pack mesh. Similar computational pa- 
rameters were successfully used in previous studies of Cu and 
Fe-doped lead titanate . 18 ' 21 ' 23 For several configurations we 
also carried out calculations using a 4 x 4 x 4 mesh and found 
negligible differences on the order of 0.05 eV and below. The 
computations were performed at the theoretical lattice con- 



stant of ao = 3.866 A and the theoretical value for the axial 
ratio of c/a = 1.05, both of which are in reasonable agree- 
ment with experiment (ao = 3.905 A, c/a — 1.064 at room 
temperature, Refs. [37] and 38). The calculated band gap of 
1.47 eV is considerably smaller than the experimental value, 
but consistent with the well known band gap error of DFT. As 
argued in Ref. |39| the band gap error is, however, expected to 
have a minor effect in the context of migration barrier calcula- 
tions. Migration paths and barriers were determined using the 
climbing image nudged elastic band method 40 ' 41 and configu- 
rations were optimized until the maximum force was less than 
30 meV/A. For charged defects a homogeneous background 
charge was added. 



B. Free oxygen vacancies 

As indicated in Fig.[T]there are two crystallographically dis- 
tinct oxygen vacancy sites in the tetragonal perovskite lattice 
(space group P4mm). Vacancies can be situated along the c- 
axis (Yc, Wyckoff site 16) or within the ao-plane {V a b, Wyck- 
off site 2c) i 18,4243 Thus, three different migration paths are 
possible between nearest neighbor sites: pure in-plane migra- 
tion (V a b —> V a b), out-of-plane migration along the positive 
direction of the tetragonal axis (V a b — > V c [001]), and out-of- 
plane migration along the negative direction of the tetragonal 
axis (V ab -> V c [001]). 

The calculated barriers are compiled in Fig. |2ja) and Ta- 
ble U In general, the barriers found to be charge state depen- 
dent, which is in line with calculations for cubic PbTi03 44 
and cubic BaTiC^? 4 ^ For migration within the a6-plane the 
barriers decrease as electrons are removed from the defect, 
which is in accord with the calculations on cubic perovskite 
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(a) free V (q=+2) (b) Cu Ti - V (q=0) 




®Fe Ti -v cl Fe Ti -l/ cl © 
(c) Fe Ti - l/ (q=+1) (d) Fe Ti - V (q=+2) 



FIG. 2. Schematic representation of the barriers for the migra- 
tion of free (a) and complexed (b-d) oxygen vacancies in units of 
eV. Each figure shows the projection of a BOe octahedron onto the 
(100) plane. The numbers in circles indicate the indices used to dis- 
tinguish the different processes. Oxygen sites (and thus possible va- 
cancy sites) are shown as red circles while the position of the B-site 
cation (Ti, Cu, or Fe), which is situated at the center of the oxygen 
octahedron, is indicated by blue (Ti), green (Cu) and yellow (Fe) cir- 
cles. 



TABLE I. Migration energies in eV for unbound oxygen vacancies 
in tetragonal lead titanate. In cases for which the jump is asymmetric 
both the backward barrier is given in brackets. 



Transition 







+1 


+2 


Vab -> Vab 


0.98 


0.62 


0.53 


V c -> Vab 


[001] 


0.83 (0.91) 


0.94 (0.58) 


1.10 (0.54) 


Vc -> Vab 


[001] 


0.50 (0.59) 


0.58 (0.22) 


0.70 (0.13) 



structures.^"* The barriers for migration via c-type vacan- 
cies, in constrast, increase. Since the charge state q = +2 
prevails vor the oxygen vacancy almost over the entire band 
gapr^ its respective barriers were used in the construction of 
the energy landscapes used in the kinetic model. 

As detailed in the appendix, we can obtain the oxygen dif- 
fusivity (excluding the formation energy contribution) from 
our calculated barriers and compare it with experimental data. 
We find that the calculated activation energy of 0.7 eV is in 
good agreement with the experimental value of 0.87 eV4* 



C. Complexed oxygen vacancies 

The incorporation of an impurity on the B site breaks trans- 
lational symmetry along the tetragonal axis and lifts the de- 
generacy of the oxygen sites in this direction. One therefore 
obtains three distinct types of first-neighbor impurity atom- 
oxygen vacancy associates (compare Figs. Q] and |2), which 
leads to three distinct migration paths. Migration within the 
second neighbor shell of impurity atoms was not considered, 
since it has been previously shown that the binding energy 
between oxygen vacancies and Cu and Fe impurities is the 
largest in the first neighbor shell. 18 Thus, once an oxygen va- 
cancy arrives in the vicinity of an impurity via diffusion, it 
will be attracted to the impurity and eventually reside in its 
first neighbor shell. The barriers for different paths for the 
migration of oxygen vacancies in the first nearest neighbor 
shell of copper and iron impurities are shown in Fig. |2}b-d). 



D. Construction of energy landscape 

So far, we have calculated the migration barriers for free 
oxygen vacancies, which determine the elements of the rate 
matrix for = {4, 5}, and the migration barriers for 
oxygen vacancies in the first neighbor shell of an impurity 
atoms, which provide the elements with = {1,2,3}. Us- 
ing these data some parts of the energy landscape can already 
be constructed as indicated by the migration barriers shown in 
black in Fig. [3] 

To determine the barriers for the remaining combinations of 
i and j, e.g., 1 — 5 or 3 — 4, one would require noticeably larger 
supercells than the ones employed in the present work. This 
results from the long ranged Coulombic attraction between 
the impurity atom and the oxygen vacancy which leads to a 
gradual transition from a free oxygen vacancy to a vacancy 
in the first neighbor shell of an impurity ion over several lat- 
tice spacings. This complexity can in principle be captured by 
increasing the dimensionality of the probability and rate ma- 
trices, $ and K. In the present more simple description, we, 
however, consider already oxygen vacancies in the second im- 
purity neighbor shell as "free". 

For completing the rate matrix, we then assume the mi- 
gration barriers for jumps between free oxygen vacancies to 
hold for jumps to complexed vacancies as well (values marked 
in green in Fig. |3). To determine the barriers for the re- 
verse jumps, we resort to the binding energies calculated for 
Cuji — Vo and Fej; — Vq complexes calculated in Ref. [j~8| 
(values marked in blue in Fig. [3}. For Cu and Fe the binding 
energy amounts to — 2.38eV and — 1.32eV for Fermi levels 
near mid gap, which leads to the energy surfaces shown in 
Fig. [3] We have tested the sensitivity of the results of the 
kinetic model to the barriers for jumps between free and com- 
plexed vacancies, which showed the assumptions made in de- 
termining their values to be of little consequence. 
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FIG. 3. Energy surface for the migration of oxygen vacancies in (a) Cu-doped (Cu-n — Vb, q = 0) and (b) Fe-doped PbTiC>3 (Fe-n — Vo, 
q — +1). The sequence of barriers corresponds to a (hypothetical) continuous trajectory which illustrates all possible migration barriers 
(compare Figs.Q]and[2j. 



IV. RESULTS 

In this section we present results obtained using the kinetic 
model for Cu and Fe-doped lead titanate described in Sect. [HI 
To illustrate the general features of vacancy redistribution, we 
first discuss the results for Cu-doped PbTiC<3 in the absence 
and presence of electric fields in Sects. |IV A| and IIVBI re- 
spectively. The results for Fe-doped material are qualitatively 
very similar and will be summarized in Sect. IIV CI 



A. Vacancy redistribution in the absence of electric fields 

Using the energy landscape for Cu-doped PbTiC>3 shown in 
Fig. |3] as input data for the kinetic model one can obtain the 
temporal evolution of the relative concentrations of different 
types of vacancies as exemplarily shown in Fig.|4f a) for a tem- 
perature of 300 K and a dopant concentration of /m = 5%. 
In this example, the vacancies are initially statistically dis- 
tributed over all available sites. In Fig. Ufa) four distinct time 
regimes with characteristically different dynamic balance can 
be identified that are separated by the transitions marked A, 
B, and C: 

(i) The first regime (up to t < 10~ 10 s at 300 K) is associ- 
ated with the redistribution of unbound vacancies. As can be 
seen from Fig. [3] V c is energetically preferred over V a b- Fig- 
ure |4ja) shows that even for a temperature as low as 300 K 
the redistribution between these two types of vacancies, i.e. 
the installation of the (partial) equilibrium over the subset of 
unbound vacancies, takes place within fractions of a second. 

(ii) During the second stage (10" 10 s - 10° s at 300 K) 
unbound c-type vacancies dominate. Concurrently, dopants 



begin to capture vacancies. The dynamic equilibrium at 
this point is such that Cuji — V C 2 complexes dominate over 
Cuxi — V c \ and Cuxi — V a b dipoles. The former are ener- 
getically less favorable but are much more easily accessible 
since AG m (V a b — > Mb — V&) is only 0.13eV compared to 
AG m (V ab -> M B - V ab ) = 0.53 eV, AG m (V c -> M B - 
V ab ) = 0.70 eV, and AG m (V ab -> M B -V cl ) = 0.54 eV. 

(Hi) In the third stage (10° s - 10 5 s at 300 K) vacancy- 
dopant complexes take over with Cu-n — V C 2 being the dom- 
inant defect. The prevalence of Cut; — V C 2 is inherited from 
the second stage which determines the initial concentrations 
for the third stage. 

(iv) Eventually, the system reaches equilibrium (t > 10 5 s 
at 300 K), i.e. virtually all vacancies occupy the lowest energy 
site£ 

We can now discuss the ependence of vacancy migra- 
tion on temperature and dopant concentration. As indi- 
cated by the letters A, B, and C in Fig. SJa), characteristic 
times can readily identified, at which the majority defect type 
changes. As shown in Fig.|4fb) the temperature dependence of 
these characteristic times can be fit to an Arrhenius equation, 
t _1 = wq exp[— Eji/k B T], using the migration barrier be- 
tween the states involved as the activation energy. This anal- 
ysis also demonstrates that the effect of changing the dopant 
concentration is small and is only visible for the transition be- 
tween vacancy types V c and Cuji — V C 2- 

Already at temperatures > 450 K the full equilibrium is 
established within less than a second. Since during growth 
these temperatures are easily reached, the vacancy distribu- 
tion in tetragonal lead titanate should be in thermal equilib- 
rium, i.e. virtually all dopants are complexed with vacancies 
in the ground state configuration, in which the defect dipoles 
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FIG. 4. Vacancy redistribution in the absence of electric fields for 
Cu-doped lead titanate: (a) Temporal evolution of the relative con- 
centrations of different vacancy types at a temperature 300 K and 
for a Cu concentration of 5%. (b) Temperature dependence of the 
characteristic transition times marked by white circles in (a). The 
transition time for the conversion from V c to Cun — V c 2 depends on 
the Cu concentration as exemplified by the three green dashed lines 
of varying thickness. 



are aligned with the domain polarization. 



B. Vacancy redistribution in the presence of electric fields 

In the present model the perturbation introduced by an 
external electric field enters in two ways. First the barri- 
ers for vacancy jumps with components along the direction 
of the electric field are distorted ("direct effect"), AG m — > 
AG m - 8E where 5E = Ei oc Ar[ 001] qe. Here, Ar [001 ] 
denotes the displacement along the direction of the electric 
field which is positive (negative) if the displacement is par- 
allel (anti-parallel) to the electric field; q is the charge state 
of the defect, Ei oc is the local electric field, and e denotes 
the unit charge. Typical electric fields used for poling fer- 
roelectric ceramics are on the order of 2-5 kV/cm; the local 
electric field can, however, be larger than this value due to 
inhomogeneities. 30 Assuming a value of Ei oc = lOOkV/cm 
and choosing Ar^oi] = 1 A, we obtain an upper limit for SE 
of 0.01 eV. 

Obviously the "direct" effect of the electric field is rather 
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FIG. 5. Vacancy redistribution in the presence of an oscillating 
external field: Equilibration over vacancy types Mb-V c i (solid lines) 
and Mb-Vc2 (dashed lines) in the presence of an external oscillating 
field with a cycling frequency of 1 Hz for temperatures between 300 
and 360 K. 



small compared to the energy difference between different va- 
cancy types and pertains to charged vacancies only. This im- 
plies that with regard to vacancy redistribution a material in 
a constant external field will behave almost identical to the 
situation without electric fields. 

The situation does, however, change if we consider an os- 
cillating external field. The field induced reversal of the polar- 
ization has a much larger effect on the energetics of the system 
than the direct contribution since reorientation of the polariza- 
tion implies that the (average) displacement of B-site atoms is 
reversed, thus transforming Mb-V c \ into Mb-V C 2 complexes. 
In the present model, this is equivalent to exchanging rows 1 
and 2 of the migration rate matrix, Ky. One can therefore 
include the effect of an oscillating electric field by (i) periodi- 
cally modifying the barriers for out-of -plane jumps by 5E and 
(z'z) simultaneously exchanging the barriers for jumps involv- 
ing Mb-V c i or M B -V c2 . 

Figure [5] illustrates the temporal evolution of the concen- 
trations of Mb-V c \ and Mb-V C 2 vacancies in the presence of 
an oscillating electric field. The plot shows that under pro- 
longed cycling a dynamic balance between the two configu- 
rations Mb-V c i and Mb-V c i is established. The time after 
which this balance is obtained depends sensitively on tem- 
perature e.g., at room temperature it is reached after about 
10 6 s (approximately two weeks) while at 340 K it requires 
only about one hour. 

The dynamic balance between Mb-V& an d Mb- V C 2 occurs 
because the characteristic time required to reach full equilib- 
rium in the absence of external electric fields [see Fig. SJb)] 
exceeds the cycling period. In a fast switching field the restor- 
ing force of the spontaneous polarization is changing signs 
on a short time scale. As a result, a mean-field composed 
of parallel and anti-parallel polarization of the matrix is act- 
ing on the dipoles, which evantually populate both orienta- 
tions along the c-axis. This dynamic equilibrium is sensitive 
to temperature and the assumption of a static distribution of 
defect dipoles in a fast switching field should be taken with 
care. The results show that by applying a bipolar electric field 
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FIG. 7. Schematic representation of defect dipole arrangements un- 
der different conditions as deduced from the kinetic model. The large 
value arrows indicates the matrix polarization in different domains 
while the small arrows represent defect dipoles. The thick solid line 
illustrates the position of a 90° domain wall. 



pare Fig. [2j, which speeds up the transition. We can thus ex- 
pect that in the presence of an oscillating external field the 
dynamic equilibrium between Feji — 142 and Fej; — V c i is es- 
tablished faster as well, which is confirmed by explicit calcu- 
lation. Whereas for Cu-doped lead titanate our model calcu- 
lations predict the equilibrium to be installed over two weeks 
at room temperature, in Fe-doped material the same process 
should occur on the order of a day. Similarly at 340 K equili- 
bration should take only on the order of tens of minutes. 



Inverse temperature (1/K) 



FIG. 6. (a) Temporal evolution of the relative concentrations of 
different vacancy types for three different temperatures in Fe-doped 
lead titante. (b) Dependence of characteristic time scales for the tran- 
sitions indicated by white circles in (a) on temperature and dopant 
concentration. 

defect dipoles are redistributed and on average the clamping 
effect on domain walls is reduced. This is in line with recent 
results on deaging of doped PZTi^S 

C. Results for Fe-doped lead titanate 

Complexes of Fe with oxygen vacancies act as donors lead- 
ing to electron chemical potentials in the upper half of the 
band gap. The most stable charge state is q = +1. 18 Under 
such conditions the binding energy is — 1.32eV from which 
one can construct an energy surface shown in Fig. 0b). 

The temporal evolution of different vacancy configurations 
in the absence of an external electric field is shown in Fig. |6ja) 
which allows us to infer the temperature dependence of the 
characteristic time scales summarized in Fig. 0b). Compar- 
ing Fig. [4] and Fig. [6] we find that the results for Fe and Cu- 
doped lead titanate are very similar. This is expected since 
the first two transitions (V a b — > V c and V c — > Fe-n — V C 2, 
compare Sect. MI Db are determined by the migration barri- 
ers in the pure host. With regard to the third transition be- 
tween Mr/j — V C 2 and Mn — V c i the situation is different as 
the effective barrier in Fe-doped material is l.OOeV and thus 
slightly smaller than in Cu-doped lead titanate (1.06 eV, com- 



V. SUMMARY AND DISCUSSION 

We have parametrized a kinetic model for defect dipole for- 
mation and switching by taking data from first-principles cal- 
culations for Fe and Cu-doped PbTiC>3. We find that at tem- 
peratures > 450 K, which is well below the Curie temperature 
of 720 K, 48 the formation and alignment of defect dipoles in 
doped PbTiC<3 should occur within less than a second. 49 

Bipolar poling leads to a dynamic equilibrium between de- 
fect dipoles that are aligned parallel and anti-parallel to the lat- 
tice polarization, respectively, and thus can be seen as one ma- 
jor contribution to deaging of PZT ceramics. This is in accord 
with the experimental observation of deaging by the applica- 
tion of AC fieldsi^ 9 ^ . In Cu-doped lead titanate the dynamic 
equilibration takes about two weeks at 300 K, but is massively 
accelerated if temperature is slightly increased. This points to 
the importance of closely monitoring the sample temperature 
during testing and studying aging and deaging processes. 

For the Mb-V c \ complex, which for both Cu and Fe is 
the ground state configuration, the local polarization is par- 
allel to the polarization of the surrounding matrix. 18 Since in 
contrast Mb-V C 2 defects are aligned anti-parallel to the lat- 
tice polarization, an increase in their concentration causes an 
overall loss of switchable polarization. This direct contribu- 
tion should scale linearly with the number of impurity atoms 
in the sample but due to the small magnitude of the defect 
dipole moment will amount to a rather small contribution on 
the macroscopic scale. Defect dipoles, however, also interact 
with domain walls and can affect their mobility. In lead ti- 
tanate and tetragonal PZT one typically observes 90° domain 
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wall configurations, which is schematically indicated in Fig. 
17] It has been shown by first-principles calculations that the 
head-to-tail domain wall configuration shown e.g., in Fig. [7] 
is energetically more stable than head-to-head or tail-to-tail 
configurations. 50 In the pristine material after cooling (middle 
panel of Fig. |7]i all defect dipoles are aligned with the lattice 
polarization and thus follow the head-to-tail pattern. This situ- 
ation changes significantly after cyclic loading (right panel of 
Fig. |7]) since now half of the defect dipoles oppose the lattice 
polarization and thus create local high-energy head-to-head 
and tail-to-tail configurations. 

Recent simulations of domain wall motion using an em- 
pirical force field 16 have provided impressive evidence that 
domain wall motion proceeds via a nucleation-and-growth 
process. It will be subject of future work to determine the 
role of defect dipoles quantitatively but already at the present 
stage one can imagine that defect dipoles that are aligned anti- 
parallel to the lattice polarization in the growing domain will 
locally impede both nucleation and growth while the oppo- 
site can be said for defect dipoles that are aligned parallel to 
the lattice polarization. One can expect that even though the 
direct contribution of defect dipoles to the macroscopic polar- 
ization is small they can have a significant indirect impact by 
pinning domain walls and reducing their mobility. It should 
be stressed that the fact that domain motion occurs via nucle- 
ation and growth is crucial in this context since it implies that 
domain wall motion occurs locally and can thus be strongly 
influenced by localized defect dipoles. 



switching kinetics of defect dipoles in ferroelectrics, which 
is relevant for understanding aging and deaging mechanisms. 
The rate equation approach can be adapated straightforwardly 
to describe more complex geometries and systems. This could 
be used to model the lattice geometries of e.g., BiFeC>3 or 
LaMn0 3 . 

Kinetic models similar to the one discussed in the present 
paper can also be used to interpret experimental measure- 
ments. To this end, one could employ probes which are sensi- 
tive to the orientation of the defect dipoles (e.g., electron spin 
resonance 21 ' 23 ) and measure the intensity of the signal before, 
during and after cycling or heat treatments. 
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Appendix: Oxygen Diffusivity 

We can employ the calculated migration barriers to derive 
the diffusivity of oxygen vacancies (see e.g., Ref. 1391). The 
rate at which a vacancy jumps along a given path i is 



VI. CONCLUSIONS 



i/ exp(-/3AGi) 



(A.l) 



In the present work we have derived a kinetic model that al- 
lows us to study the temporal evolution of the concentration of 
different types of vacancies both in the absence and presence 
of electric fields. The most important input parameter is the 
energy landscape for oxygen vacancy migration. Using pa- 
rameters for Cu and Fe-doped PbTiC»3 obtained from density 
functional theory calculations, we found that the equilibration 
of the vacancy distribution occurs readily at temperatures con- 
siderably below the Curie temperature. As a result in the as- 
synthesized material virtually all impurity atoms are associ- 
ated with vacancies forming Mb-V c \ complexes. The com- 
plete realignment of vacancy-metal impurity dipoles parallel 
to the spontaneous polarization occurs on time scales of hours 
to days at room temperature, but is massively accelerated if 
temperature is slightly increased. This provides evidence for 
the fact that aging due to defect dipoles occurs instantaneously 
in PbTiC<3-based ferroelectrics. 

In the presence of an oscillating electric field a dynamic 
balance between Mb-V c % and Mb-V C 2 is established. Pro- 
longed cycling therefore leads to the accumulation of defect 
dipoles that oppose the polarization of the encompassing do- 
main. While these defect dipoles directly reduce the switch- 
able polarization, more importantly they can impede domain 
wall motion, which has been recently shown to proceed via 
nucleation and growth*^ 

The present results provide valuable insights into the 



where v§ is the attempt frequency, AG; is the barrier which 
has to be surpassed along path i, and j3 — l/ksT. Summing 
over all paths and including the jump lengths A; as well as the 
path multiplicities (i the defect diffusivity is then given by 



D d = 



(A.2) 



There are four symmetrically equivalent migration processes 
(C = 4) within the afr-plane (V a b — > V a b), for which the 
jump lengths projected onto the a6-plane and the c-axis are 
A^ = ao and A|| = 0, respectively. With regard to of out-of- 
plane migration (V a b —> V c ) there are again four possibilities 
each for jumps with components along [001] and [001], re- 
spectively, associated with displacements Aj_ = ao/ v2 and 
A|| = ±-|co- Inserting these parameters into Eq. (I A. 2) yields 



D i 



D \\ = o c o^o 



S p(-^AGLn)+-p(-^Lf ] c ) 

f 2 exp (-(3AG ab - a b) exp (-0AG f a 
exp(- /3 AGLn)+ex P (^AGLf I c 



orm 
b-c 



where the very last term takes into account the equilibrium 
occupancy of afe-sites with respect to c-site vacancies and 
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AG a ^™ is the difference between the formation free energies 
of c and a6-site vacancies. Finally, the isotropic diffusivity is 
given as the trace of the diffusivity tensor which for tetragonal 
symmetry yields D = 2D± + Dn. 

According to Table Q] the process V c —> V a b[001] has the 
lowest barrier for all charge states (also compare the barriers 
between the four leftmost minima in Fig. [3j- Therefore, the 
isotropic defect diffusivity is approximately 

D a (2a 2 Q + ±cg) v a cxp (-/3AG^ ! C ) . (A.3) 

Using the migration barriers for charge state q = +2, approx- 
imating the attempt frequency by the lowest optical mode at 



r, h*o « 2 THz^ and using the experimental lattice constants 
one obtains 

D a 7.8 x Kr 3 exp(-0.7eV/£; B T)cm 2 /s. 

The activation barrier in this expression is in reasonable agree- 
ment with recent diffusion measurements on lead titanate- 
zirconate (PZT) alloys^, in which the migration barrier for 
oxygen vacancy migration was found to be 0.87 eV. 

It should be noted that the experiments in Ref.[46| were car- 
ried out at low temperatures on samples that contained an ex- 
trinsic concentration of oxygen vacancies. As a result, the 
calculated and measured pre-factors cannot be directly com- 
pared since the latter contains the (unknown) concentration of 
extrinsic vacancies in the samples. 
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